PROJECTO FINAL CONVECIÓN SOBRE PLACAS SOBRE PLACAS¶
Por:
- Luis Alejandro Rodríguez, Cod. 202321287
- Daniel Suaza,
- Santiago Ortegón,
Librerías¶
In [1]:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from IPython.display import HTML
from mpl_toolkits.mplot3d import Axes3D # necesario para 3D en algunos entornos
Organización¶
El notebook esta organizado por puntos, algunos puntos estan subdivididos de forma adicional a la recomendada en el documento del proyecto pues esto hacía más facil entender como funciona el código y hacerlo más legible.
1a¶
Primer intento modelar conducción¶
In [2]:
Tp = 1.0
T_inf = 0.0
alpha = 1.0
dy = 1.0
dt = 0.1
Ny = 40
Nt = 40
lam = alpha * dt / dy**2
print("lambda =", lam)
y = np.arange(Ny + 1) * dy
t = np.arange(Nt + 1) * dt
T = np.zeros((Nt + 1, Ny + 1))
T[0, :] = T_inf
T[0, 0] = Tp
for k in range(Nt):
Th = T[k, :].copy()
TE = Th.copy()
TE[0] = Tp
TE[-1] = T_inf
for n in range(1, Ny):
TE[n] = Th[n] + lam * (Th[n+1] - 2.0*Th[n] + Th[n-1])
T[k+1, :] = TE
Nx = 10
dx = 0.5
x_positions = np.arange(Nx) * dx
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
lines = []
for xi in x_positions:
X = np.full_like(y, xi, dtype=float)
line, = ax.plot(X, y, T[0, :], lw=2)
lines.append(line)
ax.set_xlabel("x (a lo largo de la placa)")
ax.set_ylabel("y (normal a la placa)")
ax.set_zlabel("T")
ax.set_xlim(x_positions[0], x_positions[-1])
ax.set_ylim(y[0], y[-1])
ax.set_zlim(T_inf, Tp * 1.05)
title = ax.set_title(f"t = {t[0]:.1f}")
def update(frame):
for i, xi in enumerate(x_positions):
X = np.full_like(y, xi, dtype=float)
lines[i].set_data(X, y)
lines[i].set_3d_properties(T[frame, :])
title.set_text(f"t = {t[frame]:.1f}")
return lines + [title]
anim = FuncAnimation(fig, update, frames=Nt+1, interval=200, blit=False)
plt.close(fig)
HTML(anim.to_jshtml())
lambda = 0.1
Out[2]:
Intentar cambiar orientación del modelo¶
In [3]:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from IPython.display import HTML
from mpl_toolkits.mplot3d import Axes3D # a veces necesario para activar 3D
Tp = 1.0
T_inf = 0.0
alpha = 1.0
dy = 1.0
dt = 0.1
Ny = 40
Nt = 40
lam = alpha * dt / dy**2
print("lambda =", lam)
y = np.arange(Ny + 1) * dy
t = np.arange(Nt + 1) * dt
T = np.zeros((Nt + 1, Ny + 1))
T[0, :] = T_inf
T[0, 0] = Tp
for k in range(Nt):
Th = T[k, :].copy()
TE = Th.copy()
TE[0] = Tp
TE[-1] = T_inf
for n in range(1, Ny):
TE[n] = Th[n] + lam * (Th[n+1] - 2.0*Th[n] + Th[n-1])
T[k+1, :] = TE
Nx = 10
dx = 0.5
x_positions = np.arange(Nx) * dx
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
lines = []
for xi in x_positions:
X = np.full_like(y, xi, dtype=float) # eje x (horizontal)
Y = T[0, :] # eje y (vertical) = T
Z = y # eje z (profundidad) = distancia normal
line, = ax.plot(X, Y, Z, lw=2)
lines.append(line)
ax.set_xlabel("x (a lo largo de la placa)")
ax.set_ylabel("T")
ax.set_zlabel("y (normal a la placa)")
ax.set_xlim(x_positions[0], x_positions[-1])
ax.set_ylim(T_inf, Tp * 1.05)
ax.set_zlim(y[0], y[-1])
ax.view_init(elev=25, azim=-60) # solo para que se vea más claro
title = ax.set_title(f"t = {t[0]:.1f}")
def update(frame):
for i, xi in enumerate(x_positions):
X = np.full_like(y, xi, dtype=float)
Y = T[frame, :] # T(y,t) en el eje vertical
Z = y
lines[i].set_data(X, Y) # x, T
lines[i].set_3d_properties(Z) # y normal
title.set_text(f"t = {t[frame]:.1f}")
return lines + [title]
anim = FuncAnimation(fig, update, frames=Nt+1, interval=200, blit=False)
plt.close(fig)
HTML(anim.to_jshtml())
lambda = 0.1
Out[3]:
Sub-dividir en diferentes gráficas¶
In [4]:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from IPython.display import HTML
Tp = 1.0
T_inf = 0.0
alpha = 1.0
dy = 1.0
dt = 0.1
Ny = 40
Nt = 40
lam = alpha * dt / dy**2
print("lambda =", lam)
y = np.arange(Ny + 1) * dy
t = np.arange(Nt + 1) * dt
T = np.zeros((Nt + 1, Ny + 1))
T[0, :] = T_inf
T[0, 0] = Tp
for k in range(Nt):
Th = T[k, :].copy()
TE = Th.copy()
TE[0] = Tp
TE[-1] = T_inf
for n in range(1, Ny):
TE[n] = Th[n] + lam * (Th[n+1] - 2.0*Th[n] + Th[n-1])
T[k+1, :] = TE
Nx_panels = 5
x_positions = np.linspace(0.0, 1.0, Nx_panels)
fig, axes = plt.subplots(1, Nx_panels, figsize=(4 * Nx_panels, 4), sharey=True)
lines = []
for i, ax in enumerate(axes):
line, = ax.plot(y, T[0, :], lw=2)
lines.append(line)
ax.set_xlabel("y (normal a la placa)")
ax.set_title(f"x = {x_positions[i]:.2f}")
axes[0].set_ylabel("T(y,t)")
fig.suptitle(f"t = {t[0]:.1f}")
def update(frame):
for line in lines:
line.set_data(y, T[frame, :])
fig.suptitle(f"t = {t[frame]:.1f}")
return lines
anim = FuncAnimation(fig, update, frames=Nt+1, interval=200, blit=False)
plt.close(fig)
HTML(anim.to_jshtml())
lambda = 0.1
Out[4]:
Integración final y resultado final¶
In [5]:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from IPython.display import HTML
Tp = 1.0
T_inf = 0.0
alpha = 1.0
dy = 1.0
dt = 0.1
Ny = 40
Nt = 40
lam = alpha * dt / dy**2
print("lambda =", lam)
y = np.arange(Ny + 1) * dy
t = np.arange(Nt + 1) * dt
T = np.zeros((Nt + 1, Ny + 1))
T[0, :] = T_inf
T[0, 0] = Tp
for k in range(Nt):
Th = T[k, :].copy()
TE = Th.copy()
TE[0] = Tp
TE[-1] = T_inf
for n in range(1, Ny):
TE[n] = Th[n] + lam * (Th[n+1] - 2.0*Th[n] + Th[n-1])
T[k+1, :] = TE
Nx_profiles = 10
offsets = np.linspace(0.0, 3.0, Nx_profiles)
fig, ax = plt.subplots(figsize=(7, 5))
lines = []
for off in offsets:
line, = ax.plot(off + T[0, :], y, lw=2)
lines.append(line)
ax.set_xlabel("x (perfiles repetidos a lo largo de la placa)")
ax.set_ylabel("y (normal a la placa)")
ax.set_xlim(offsets[0] - 0.2, offsets[-1] + Tp + 0.2)
ax.set_ylim(y[0], y[-1])
title = ax.set_title(f"t = {t[0]:.1f}")
def update(frame):
for i, off in enumerate(offsets):
lines[i].set_data(off + T[frame, :], y)
title.set_text(f"t = {t[frame]:.1f}")
return lines + [title]
anim = FuncAnimation(fig, update, frames=Nt+1, interval=200, blit=False)
plt.close(fig)
HTML(anim.to_jshtml())
lambda = 0.1
Out[5]:
1b¶
In [6]:
import numpy as np
import matplotlib.pyplot as plt
Tp = 1.0
T_inf = 0.0
alpha = 1.0
dy = 1.0
dt = 0.1
Ny = 40
Nt = 40
lam = alpha * dt / dy**2
print("lambda =", lam)
y = np.arange(Ny + 1) * dy
t = np.arange(Nt + 1) * dt
T = np.zeros((Nt + 1, Ny + 1))
T[0, :] = T_inf
T[0, 0] = Tp
for k in range(Nt):
Th = T[k, :].copy()
TE = Th.copy()
TE[0] = Tp
TE[-1] = T_inf
for n in range(1, Ny):
TE[n] = Th[n] + lam * (Th[n+1] - 2.0*Th[n] + Th[n-1])
T[k+1, :] = TE
k_f = 1.0
dTdy_0 = np.zeros(Nt + 1)
h = np.zeros(Nt + 1)
for k in range(Nt + 1):
dTdy_0[k] = (-1.5*T[k, 0] + 2.0*T[k, 1] - 0.5*T[k, 2]) / dy
h[k] = -k_f * dTdy_0[k] / (Tp - T_inf)
plt.figure()
plt.plot(t, h, marker="o")
plt.xlabel("tiempo t")
plt.ylabel("h(t)")
plt.title("Coeficiente de convección h vs tiempo (difusión pura)")
plt.grid(True)
plt.show()
lambda = 0.1
2a y 2b (2a esta respondido en documento informe)¶
In [7]:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from IPython.display import HTML
Tp = 1.0
T_inf = 0.0
alpha = 1.0
dy = 1.0
dt = 0.1
Ny = 40 # nodos en y (0..40)
Nt = 50 # pasos de tiempo
Nx = 40 # "puntos" en x
U_inf = 1.0
dx = U_inf * dt # para ser consistente con T(x+U dt)
lam = alpha * dt / dy**2
print("lambda =", lam)
y = np.arange(Ny + 1) * dy
t = np.arange(Nt + 1) * dt
x = np.arange(Nx) * dx
# T[k, i, n] = T en tiempo k, posición x_i, nodo y_n
T = np.zeros((Nt + 1, Nx, Ny + 1))
# Condición inicial: fluido frío, placa caliente en todo x
T[0, :, :] = T_inf
T[0, :, 0] = Tp
for k in range(Nt):
# -------- Paso 1: CONDUCCIÓN en y para cada x_i --------
T_cond = np.zeros_like(T[k])
for i in range(Nx):
Th = T[k, i, :].copy()
TE = Th.copy()
TE[0] = Tp # placa
TE[-1] = T_inf # medio lejano
for n in range(1, Ny):
TE[n] = Th[n] + lam * (Th[n+1] - 2.0*Th[n] + Th[n-1])
T_cond[i, :] = TE
# -------- Paso 2: TRANSPORTE en x --------
# Queremos:
# - Columna i=0: siempre fluido frío arriba y placa a Tp (estática)
# - Columna i>=1: se llena con lo que había en i-1 después de la conducción
T_next = np.zeros_like(T_cond)
# i = 0: condición de entrada (fluido nuevo)
T_next[0, :] = T_inf
T_next[0, 0] = Tp
# i >= 1: transportar perfiles del fluido (y>0) desde la columna anterior
for i in range(1, Nx):
# copiamos solo el fluido (n>=1) desde i-1
T_next[i, 1:] = T_cond[i-1, 1:]
# placa siempre a Tp
T_next[i, 0] = Tp
T[k+1, :, :] = T_next
# -------- ANIMACIÓN: curvas lado a lado sobre el eje x --------
offsets = x # usamos x como desplazamiento horizontal
fig, ax = plt.subplots(figsize=(10, 5))
lines = []
for i in range(Nx):
# al inicio k=0
line, = ax.plot(offsets[i] + T[0, i, :], y, lw=0.7)
lines.append(line)
ax.set_xlabel("x (a lo largo de la placa)")
ax.set_ylabel("y (normal a la placa)")
ax.set_xlim(offsets[0] - 0.2, offsets[-1] + Tp + 0.5)
ax.set_ylim(y[0], y[-1])
title = ax.set_title(f"t = {t[0]:.1f}")
def update(frame):
for i in range(Nx):
lines[i].set_data(offsets[i] + T[frame, i, :], y)
title.set_text(f"t = {t[frame]:.1f}")
return lines + [title]
anim = FuncAnimation(fig, update, frames=Nt+1, interval=120, blit=False)
plt.close(fig)
HTML(anim.to_jshtml())
lambda = 0.1
Out[7]:
3a¶
In [8]:
import numpy as np
import matplotlib.pyplot as plt
# ---------------- PARÁMETROS ----------------
Tp = 1.0
T_inf = 0.0
alpha = 1.0
dy = 1.0
dt = 0.1
Ny = 40 # nodos en y (0..40)
Nt = 50 # pasos de tiempo (k = 0..50)
Nx = 40 # nodos en x
U_inf = 1.0
dx = U_inf * dt # Δx = U_inf * Δt
lam = alpha * dt / dy**2
print("lambda =", lam)
y = np.arange(Ny + 1) * dy
t = np.arange(Nt + 1) * dt
x = np.arange(Nx) * dx
# T[k, i, n] = T en tiempo k, posición x_i, nodo y_n
T = np.zeros((Nt + 1, Nx, Ny + 1))
# Condición inicial: fluido frío, placa caliente en todo x
T[0, :, :] = T_inf
T[0, :, 0] = Tp
# ---------------- CONDUCCIÓN + TRANSPORTE ----------------
for k in range(Nt):
# Paso 1: CONDUCCIÓN en y para cada x_i
T_cond = np.zeros_like(T[k])
for i in range(Nx):
Th = T[k, i, :].copy()
TE = Th.copy()
TE[0] = Tp # placa
TE[-1] = T_inf # medio lejano
for n in range(1, Ny):
TE[n] = Th[n] + lam * (Th[n+1] - 2.0*Th[n] + Th[n-1])
T_cond[i, :] = TE
# Paso 2: TRANSPORTE en x
# i = 0: fluido nuevo (frío) + placa caliente
T_next = np.zeros_like(T_cond)
T_next[0, :] = T_inf
T_next[0, 0] = Tp
# i >= 1: copiar fluido (n>=1) desde la columna anterior, placa fija a Tp
for i in range(1, Nx):
T_next[i, 1:] = T_cond[i-1, 1:]
T_next[i, 0] = Tp
T[k+1, :, :] = T_next
# ---------------- PARTE 3a: h(x) EN k = 50 ----------------
time_index = 50 # k = 50
k_f = 1.0 # k_f = 1 → h = - dT/dy / (Tp - T_inf)
dTdy_0 = np.zeros(Nx)
h_x = np.zeros(Nx)
for i in range(Nx):
T0 = T[time_index, i, 0]
T1 = T[time_index, i, 1]
T2 = T[time_index, i, 2]
dTdy_0[i] = (-1.5 * T0 + 2.0 * T1 - 0.5 * T2) / dy # ∂T/∂y|_{y=0}
h_x[i] = -k_f * dTdy_0[i] / (Tp - T_inf) # h(x_i)
plt.figure(figsize=(7, 4))
plt.plot(x, h_x, marker="o")
plt.xlabel("x")
plt.ylabel("h(x) en k = 50")
plt.title("h(x) = - (∂T/∂y|_{y=0}) / (Tp - T_inf) para k = 50")
plt.grid(True)
plt.show()
lambda = 0.1
3b¶
In [9]:
import numpy as np
import matplotlib.pyplot as plt
# Usamos solo x > 0
mask = x > 0.0
x_pos = x[mask]
h_pos = h_x[mask]
inv_sqrt_x = 1.0 / np.sqrt(x_pos)
plt.figure(figsize=(7, 4))
plt.plot(inv_sqrt_x, h_pos, "-o") # línea + puntos unidos
plt.xlabel("1 / sqrt(x)")
plt.ylabel("h")
plt.title("h vs 1/√x (k = 50)")
plt.grid(True)
plt.show()
# (opcional) Ajuste lineal
coef = np.polyfit(inv_sqrt_x, h_pos, 1)
a, b = coef
print("Ajuste h ≈ a * (1/√x) + b")
print("a =", a)
print("b =", b)
Ajuste h ≈ a * (1/√x) + b a = 0.4491322922017203 b = 0.11117425956630435
3c Note que es la misma¶
In [10]:
import numpy as np
import matplotlib.pyplot as plt
# ---------------- PARÁMETROS COMUNES ----------------
Tp = 1.0
T_inf = 0.0
alpha = 1.0
dy = 1.0
dt = 0.1
Ny = 40
Nt = 50 # k = 0..50
Nx = 40 # nodos en x
U_inf = 1.0
dx = U_inf * dt # Δx = U_inf * Δt
lam = alpha * dt / dy**2
print("lambda =", lam)
y = np.arange(Ny + 1) * dy
t = np.arange(Nt + 1) * dt
x = np.arange(Nx) * dx
# =====================================================
# 1) SOLO CONDUCCIÓN EN Y (SEMIESPACIO) → h_cond(t)
# =====================================================
T_cond1D = np.zeros((Nt + 1, Ny + 1))
T_cond1D[0, :] = T_inf
T_cond1D[0, 0] = Tp
for k in range(Nt):
Th = T_cond1D[k, :].copy()
TE = Th.copy()
TE[0] = Tp
TE[-1] = T_inf
for n in range(1, Ny):
TE[n] = Th[n] + lam * (Th[n+1] - 2.0*Th[n] + Th[n-1])
T_cond1D[k+1, :] = TE
k_f = 1.0
h_cond_t = np.zeros(Nt + 1)
for k in range(Nt + 1):
T0 = T_cond1D[k, 0]
T1 = T_cond1D[k, 1]
T2 = T_cond1D[k, 2]
dTdy0 = (-1.5*T0 + 2.0*T1 - 0.5*T2) / dy
h_cond_t[k] = -k_f * dTdy0 / (Tp - T_inf)
# =====================================================
# 2) CONDUCCIÓN + TRANSPORTE → h_conv(x) en k = 50
# =====================================================
T = np.zeros((Nt + 1, Nx, Ny + 1))
T[0, :, :] = T_inf
T[0, :, 0] = Tp
for k in range(Nt):
# Paso 1: conducción en y
T_cond = np.zeros_like(T[k])
for i in range(Nx):
Th = T[k, i, :].copy()
TE = Th.copy()
TE[0] = Tp
TE[-1] = T_inf
for n in range(1, Ny):
TE[n] = Th[n] + lam * (Th[n+1] - 2.0*Th[n] + Th[n-1])
T_cond[i, :] = TE
# Paso 2: transporte en x
T_next = np.zeros_like(T_cond)
T_next[0, :] = T_inf
T_next[0, 0] = Tp
for i in range(1, Nx):
T_next[i, 1:] = T_cond[i-1, 1:]
T_next[i, 0] = Tp
T[k+1, :, :] = T_next
time_index = 50
h_conv_x = np.zeros(Nx)
for i in range(Nx):
T0 = T[time_index, i, 0]
T1 = T[time_index, i, 1]
T2 = T[time_index, i, 2]
dTdy0 = (-1.5*T0 + 2.0*T1 - 0.5*T2) / dy
h_conv_x[i] = -k_f * dTdy0 / (Tp - T_inf)
# =====================================================
# 3) COMPARACIÓN: h_cond(t=x/U_inf) vs h_conv(x)
# =====================================================
# x = i*dx, t_equiv = x/U_inf = i*dt → índice de tiempo ~ i
# usamos i = 1..Nx-1 para evitar x=0
i_min = 1
i_max = min(Nx-1, Nt) # por si Nx > Nt
x_common = x[i_min:i_max+1]
h_cond_on_x = h_cond_t[i_min:i_max+1] # h_cond(t = i*dt)
h_conv_on_x = h_conv_x[i_min:i_max+1]
plt.figure(figsize=(8, 4))
plt.plot(x_common, h_cond_on_x, "-o", label="Solo conducción (h(t), t=x/U_inf)")
plt.plot(x_common, h_conv_on_x, "-s", label="Conducción + transporte (h(x), k=50)")
plt.xlabel("x")
plt.ylabel("h")
plt.title("Comparación de h: solo conducción vs conducción + transporte")
plt.grid(True)
plt.legend()
plt.show()
lambda = 0.1